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Abstract:  We  present  a  Monte  Carlo  solution  to  the  distributed  data  fusion  problem  and 
apply  it  to  distributed  particle  filtering.  The  consensus-based  fusion  algorithm  is  iterative  and  it 
involves  the  exchange  and  fusion  of  empirical  posterior  densities  between  neighbouring  agents. 
As  the  fusion  method  is  Monte  Carlo  based  it  is  naturally  applicable  to  distributed  particle 
filtering.  Furthermore,  the  fusion  method  is  applicable  to  a  large  class  of  networks  including 
networks  with  cycles  and  dynamic  topologies.  We  demonstrate  both  distributed  fusion  and 
distributed  particle  filtering  by  simulating  the  algorithms  on  randomly  generated  graphs. 


1.  INTRODUCTION 

As  modern  sensor  systems  move  towards  distributed  archi¬ 
tectures  it  is  important  to  consider  the  design  of  practical 
algorithms  for  data  estimation  and  fusion.  Such  algorithms 
should  be  robust  and  promote  network  scalability. 

In  this  paper  we  develop  a  distributed  method  for  data 
fusion  over  networks  and  apply  the  algorithm  to  dis¬ 
tributed  particle  filtering.  Our  approach  to  this  problem 
is  consensus-based  and  involves  the  local  exchange  of 
posterior  densities.  The  posteriors  are  exchanged  via  sets 
of  unweighted  particles.  The  algorithm  converges  to  the 
optimal  (centralised)  solution  if  the  likelihood  functions 
are  conditionally  independent.  If  they  are  dependent  then 
our  algorithm  can  be  modified  to  return  a  posterior  that 
is  guaranteed  to  be  consistent  [Bailey  et  al.  (2012)]. 

There  are  many  attractive  aspects  to  our  approach,  (i) 
the  method  is  distributed:  the  agents  share  and  compute 
with  local  knowledge,  (ii)  the  fusion  scheme  is  robust  to 
networks  containing  cycles  and  networks  with  time  varying 
topologies,  (iii)  the  algorithm  has  a  guaranteed  speed  of 
convergence,  (iv)  the  resulting  posterior  distribution  is 
guaranteed  to  be  consistent  as  long  as  the  dependence 
between  the  agents  is  understood. 

The  algorithm  we  present  is  falls  under  the  category  of 
consensus-based  distributed  particle  filtering.  An  alterna¬ 
tive  to  consensus  is  message  passing  (or  routing)  based 
information  sharing.  The  following  is  a  short  review  of 
some  of  the  literature  in  this  space. 

The  distributed  fusion  method  of  Hlinka  et  al.  (2012)  (see 
also  Hlinka  et  al.  (2011)  and  Hlinka  et  al.  (2010))  relies  on 
the  individual  agent’s  likelihood  functions  belonging  to  the 
exponential  family  of  distributions.  They  apply  a  consen¬ 
sus  scheme  to  approximations  of  the  exponential  family 
and  arrive  at  a  joint  likelihood  function.  The  consensus 
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method  of  Gu  (2007)  involves  sharing  the  parameters  of 
Gaussian  mixture  models.  The  approaches  of  Sheng  et  al. 
(2005)  and  Hlinka  et  al.  (2009)  also  transmit  the  param¬ 
eters  of  a  Gaussian  mixture  model.  However  their  fusion 
schemes  are  not  based  upon  consensus,  but  rather  message 
passing.  Gu  et  al.  (2008),  another  consensus-based  fusion 
scheme,  simplifies  the  parameter  exchange  by  transmitting 
only  the  mean  and  covariance.  Consequently  this  supports 
Gaussian  dynamic  systems,  however  they  propose  its  use 
for  systems  with  higher  order  moments  by  applying  the 
unscented  transform.  Mohammadi  and  Asif  (2011a)  also 
take  advantage  of  the  unscented  transform  and  propose  its 
use  to  incorporate  distributed  Kalman  filters  into  networks 
of  distributed  particle  filters. 

Mohammadi  and  Asif  (2011b)  propose  another  consensus 
scheme  where  each  node  runs  two  particle  filters:  a  local 
particle  filter,  and  a  fusion  filter,  where  the  fusion  filter 
computes  the  global  distribution.  Ustebay  et  al.  (2011)  use 
a  selective  gossip  procedure  to  arrive  at  a  consensus  about 
the  likelihoods  of  particles.  The  selective  gossip  procedure 
shares  particles  based  upon  weights.  This  focuses  commu¬ 
nication  on  the  particles  that  contain  the  most  informa¬ 
tion.  Oreshkin  and  Coates  (2010)  also  use  a  gossip  con¬ 
sensus  approach  to  share  Gaussian  approximations  of  the 
posterior.  Recent  work  by  Lee  and  West  (2013)  proposes  a 
Markov  Chain  Distributed  Particle  Filter  (see  additionally 
Lee  and  West  (2009)  and  Lee  and  West  (2010))  in  which 
neighbours  exchange  particles  and  weights  using  a  Markov 
chain  random  walk. 

Further  approaches  to  distributed  particle  filtering  are 
outlined  in  Bashi  et  al.  (2003),  Sheng  and  Hu  (2005),  Bolic 
et  al.  (2005),  Mohammadi  and  Asif  (2009),  Farahmand 
et  al.  (2011),  and  Savic  et  al.  (2012).  In  addition,  the 
review  by  Hlinka  et  al.  (2013)  is  an  excellent  resource  for 
comparing  methods  of  distributed  particle  filtering. 

This  paper  is  organised  as  follows.  Section  2  briefly  dis¬ 
cuses  the  distributed  data  fusion  problem,  then  details  the 
proposed  Monte  Carlo  based  fusion  algorithm.  The  algo¬ 
rithm  is  supported  by  two  simulations.  Section  3  considers 
the  application  of  the  previously  derived  distributed  Monte 
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Carlo  fusion  algorithm  to  the  problem  of  distributed  par¬ 
ticle  filtering.  In  Section  4  we  offer  conclusions. 

2.  DISTRIBUTED  MONTE  CARLO  DATA  FUSION 


Taking  the  exponential  gives 


Q  oc  p(x)g{yi\x)1/n 


(5) 


Consider  a  group  of  agents  indexed  in  V  =  {1, . . . ,  n}  and  a 
set  of  possible  time- varying  undirected  links  £(k)  C  V  x  V 
defining  a  network  graph  Q(V1£(k)).  The  neighbor  set  at 
agent  i  is  denoted  by  Afi(k)  =  {j  €  V  :  ( i,j )  €  £(k)}  and 
j  €  Afi(k)  <t=>  *  €  A fj{k)  for  undirected  topologies.  We  often 
drop  the  network’s  dependence  on  k  for  brevity. 


Each  agent  constructs  an  initial  local  posterior  of  the  form 
g(yi\x)p{x) 


P(x\yi )  = 


I  g{Vi\x)p(x)dx 


g(yi\x)p(x) 


(1) 


from  given  measurements  yi  G  trasi  and  where  g(y%\x)  is 
the  likelihood  function  at  agent  i  conditioned  on  some  un¬ 
derlying  event  x  G  Rmx .  Here,  p{x )  is  the  prior  information 
common  to  all  agents.  The  goal  of  distributed  data  fusion 
in  this  case  is  to  compute 


P{x\{yi}iev)  ocp(x)  Y[g(yi\x)  (2) 

i£V 


locally  at  each  agent  i  under  the  constraint  that  agent 
i  can  only  share  p(x\yi)  with  its  neighbours  in  A/)(f).  In 
other  words,  each  agent  is  constrained  to  computations 
involving  local  posteriors,  p(x\yi)  and  p{x\yj)  where  j  G 
A fi-  Obviously,  an  iterative  procedure  is  required  to  reach 
p(x\{yi}i£v)  at  each  agent  in  an  incomplete  network. 


The  following  theorem  is  a  slight  modification  of  the  main 
result  in  Olfati-Saber  et  al.  (2006). 


Theorem  1.  Consider  a  network  Q(t)  as  described  above 
where  each  agent  exchanges  nlk  with  irlk=0  =  p{x\yi).  Here 
nk  and  7 rJk  are  not  conditionally  independent.  Suppose 
p(x)  >  0  and  g(yi\x)  >  0  for  all  i  G  V.  Then  the  following 
statements  hold: 


which  completes  the  proof.  □ 

Note  that  despite  the  fact  each  agent  is  sharing  its  local 
posterior  the  consensus-based  distributed  fusion  algorithm 
just  depicted  converges  to  p(x)  Kljev  giVjlx)1^  and  n°t 
P{x)n  Ujev  9(yj\x)1/n  orp(x)1/”  Hjev  g{yj\x)1/n.  Thus  it 
is  not  ‘double  counting’  the  common  prior  information 
p{x)  and  nor  is  it  conservative  in  this  common  prior 
information.  The  algorithm  is  actually  conservative  in 
Tljev  4Vi \x)1^n  which  may  be  important  for  consistency 
as  discussed  in  Bailey  et  al.  (2012).  In  particular,  if  yt 
and  yj  were  correlated  (i.e.  g{yi\x)  and  g(yj \x)  were  not 
conditionally  independent)  and  this  dependence  was  not 
considered  then  over-confident  results  may  be  obtained 
if  one  simply  computes  Iljev g(Vj\x)-  Similarly,  if  p(x) 
was  counted  multiple  times  then  over-confident  results  are 
clearly  obtained.  The  proposed  algorithm  converges  to  the 
log- linear  opinion  pool  [see  Abbas  (2009);  Bailey  et  al. 
(2012)]  on  the  likelihoods  multiplied  by  the  common  prior 
which  is  guaranteed  consistent.  It  converges  to  this  value  in 
a  distributed  manner  which  generalises  Bailey  et  al.  (2012). 

Corollary  2.  Suppose  that  g(jji\x)  and  g(yj\x)  are  con¬ 
ditionally  independent  for  all  i,  j  G  V.  Define  7r(.=0  oc 
p(x)g(yi\x)n.  The  proposed  consensus-based  distributed 
fusion  algorithm  converges  to  Q'  oc  p{x)Y\ieV  g(yi\x)  as 
k  — >  oo.  This  is  exactly  the  optimal  centralised  Bayesian 
result  given  p(x\yi)  oc  g(y.i\x)p(x)  where  g{yi\x)  and 
g(yj\x)  are  conditionally  independent  and  p(x)  is  the  com¬ 
mon  prior  information  among  all  agents. 

The  algorithm  can  be  rewritten  in  the  form 


i  The  agents  are  capable  of  asymptotically  reaching  a 
consensus  on  Q  oc  p(x)  Tljev  5(2/j \x)1^n- 

ii  The  consensus  algorithm  for  agreement  in  the  value  Q 
takes  the  form 

4  =  n  K-i)7-  (3) 

jeM 

where  0  <  7  <  l/max({|A/)|  :  i  G  V})  2  . 


Proof.  It  is  shown  in  Olfati-Saber  et  al.  (2006)  that 
7r 1  ~ ^  ri/evM)1^”  as  k  — >  00.  In  other  words,  log(7r(.)  — > 
(Eiev  l°s(4))/n  as  k  00. 

Note  then  in  this  case  7r(,=0  =  p(x\yi)  oc  g(yi\x)p(x)  and 


!og(7Tfc)  -> 

OC 


Ejev  ^(^o) 

n 

n\og(p(x))  Yyjev^°s(g(yj\x)) 
n  n 


(4) 


2  The  restriction  that  7  is  less  than  the  inverse  of  the  maximum 
degree  in  the  network  is  sufficient  for  each  agent  to  converge  to  Q. 
See  Olfati-Saber  et  al.  (2006)  or  Olfati-Saber  and  Murray  (2004). 


lk- 


in 

jeAfi 


(6) 


which  shows  that  when  consensus  is  reached  the  iterations 
reduce  to  7r(.  =  nk_1. 

This  algorithm  has  many  appealing  points.  Firstly,  it  is 
distributed:  the  agents  only  share  and  compute  with  local 
knowledge.  Secondly,  this  distributed  fusion  algorithm  can 
deal  with  cycles  and  time-varying  topologies  (where,  for 
example,  the  network  may  be  disconnected  for  many  time 
steps).  Thirdly,  the  algorithm  has  a  guaranteed  speed  of 
convergence  that  is  characterised  by  the  network’s  alge¬ 
braic  connectivity;  see  Olfati-Saber  et  al.  (2006).  Finally, 
as  noted  above,  the  proposed  algorithm  converges  to  the 
log-linear  opinion  pool  multiplied  by  the  common  prior  in¬ 
formation  (which  is  a  conservative,  guaranteed  consistent, 
version  of  the  optimal  Bayesian  result  regardless  of  the 
dependence  between  g(yi\x)  and  g{yjj\x)).  Moreover,  with 
a  slight  modification  to  the  initially  shared  data  we  can 
also  achieve  distributed  convergence  to  the  true  optimal 
Bayesian  fusion  result  as  noted  in  the  corollary  (and  this 
result  is  obviously  the  most  desired  when  it  is  known 
g{yi\x)  and  g(yj\x)  are  conditionally  independent). 
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2.1  A  Monte  Carlo  Implementation  of  the  Distributed 
Fusion  Algorithm 


Given  the  desirable  nature  of  the  proposed  data  fusion 
algorithm  it  remains  to  establish  a  Monte  Carlo  version 
which  allows  one  to  begin  with  an  empirical  estimate  of 
4=0  P{x)g{Vi\x)n  or  nlk=0  oc  p(x)g{yi\x)  given  by 

i  N° 

4= oO)  =  Y  s(x  ~ x 4  (?) 

s  e-i 

where  27  G  are  a  set  of  independent  and  identically 
distributed  samples  of  p(x)g(yi\x)n  or  p(x)g(yi\x)  etc. 


The  main  contribution  of  this  section  is  a  Monte  Carlo  ver¬ 
sion  of  the  algorithm  outlined  in  the  previous  subsection 
that  allows  one  to  begin  with  7fJ._0.  This  is  desirable  when 
dealing  with  complex  estimation  problems  and/or  when 
7r/_0  is  sourced  from  alternative  Monte  Carlo  estimations 
as  in  Doucet  et  al.  (2001). 


So  given  4_i)  Vj  €  A/)  U  {i}  we  want  to  compute  nlk  at 
each  agent  i  in  the  sense  that  4  should  be  an  empirical 
version  of 


4  =  4-1  LI 

jeMi 


(8) 


where  0  <  7  <  1/  max({|Al|  :  i  G  V}). 

If  the  supports  of  4-i>  Vj  G  A/)  U  {*}  were  totally 
overlapping  (i.e.  if  these  estimates  were  constructed  from 
the  same  sample  points)  then  one  could  simply  compute  7f). 
directly  via  (8).  This  is  essentially  the  situation  considered 
in  Savic  et  al.  (2012)  and  Lindberg  et  al.  (2013)  where  the 
initial  posteriors  4=0  are  sampled  at  the  same  point  for 
each  iG  V.  However,  in  the  more  likely  case  in  which  the 
supports  of  7f/_1,  Vj  G  Mi  U  {?'}  are  totally  disjoint  (with 
probability  1)  then  an  alternative  method  of  computing  an 
empirical  estimate  rtlk  of  nk  is  needed. 

Note  that  we  obviously  cannot  sample  directly  from  nk 
even  if  we  know  4_i>  Vj  G  A/j  U  {*}.  We  can  estimate 
4_  1,  Vj  GJV.U  {?'}  from  4_i>  Vj  G  A/)  U  {*},  using  for 
example  Kernel  density  estimation  as  in  Silverman  (1986), 
but  we  cannot  then  use  this  to  compute  an  estimate  of  nk 
as  no  closed  form  solution  to  the  update  in  (8)  exists  for 
any  reasonable  choice  of  the  Kernel  density  estimate. 

Instead  we  propose  the  following  empirical  approximation 
of  nk  based  on  importance  sampling;  see  Doucet  et  al. 
(2001).  Note 


7 r,.  =  7r 


=4-1^  n 

qk -1  jeMi 


‘fc-1 


‘fc-1 


(9) 


where  here  Xi  G  Rmx  are  a  set  of  Ns  independent  and 
identically  distributed  samples  of  the  so-called  importance 
function  qlk_1.  Two  matters  remain.  Firstly,  one  needs 
an  importance  function  qk_1  from  which  one  can  easily 
sample  from.  Secondly,  we  still  don’t  know  nk_1,  Vj  G  A/jU 

{?'}  and  thus  cannot  compute  4-i(4)>  Vj  G  A/j  U  {?'}. 
To  resolve  the  second  matter  we  resort  to  Kernel  density 
estimation  and  obtain 


N. 


^)=E!HSn  fi=^)  tf(*-*f)(u) 


frf  4- i(4)  \4-i(4) 

where  again  x\  G  Rmx  are  samples  of  qlk_1  and 

|  Ns  | 

44)  =  —Y~K 

Nx  h 


(12) 


where  here  27  G  Rmx  are  samples  of  4  or  in  other  words 
correspond  to  the  support  of  7f k  ( x )  and  thus  are  obviously 
known  (i.e.  by  assumption  at  k  =  0  we  know  Wk(x)).  The 
Kernel  K(-)  is  chosen  to  be  Gaussian  in  our  simulations 
but  other  choices  are  possible;  see  Silverman  (1986). 


The  importance  function  qk_1(x)  must  now  be  chosen  and 
given  the  information  available  at  the  agents  one  obvious 
choice  is 

4-i0)  =  4-44  (13) 

where  the  support  of  4-i(4  is  distributed  according  to 
qk_i(x)  =  4-1(4  and  so  sampling  from  4-i(4  is  given. 
In  this  case 


where 


Ns 

nl(x)  =  Ywl,e6(x  ~  xe) 

t=i 


=  4-1  (4)  TT  ( 4-l(4)V 
4-i(4)  \4- i(®j) ) 

tt  f 4-i(4A 

-  hk  Vi-M)) 


(14) 


(15) 

(16) 


with  wlke  =  wk  if  wk  i  an(l  xe  e  sampled  from 
4-i(4  =  4-i(4  which  in  practice  means  x\  are  exactly 
the  samples  of  4-i(4  at  the  previous  time  step.  To 
combat  possible  degeneracy  we  then,  given  the  wk  e  and  x\ 

that  define  Jfk(x),  resample  Ns  times  from  a  multinomial 
distribution  defined  by  wk  e  and  x\  to  obtain  the  final 

1  N° 

4(4  =  (17) 

s  t=i 


where  now  37  G  Rmx  corresponds  to  the  set  of  independent 
and  identically  distributed  samples  of  the  multinomial 
distribution  just  discussed. 


2.2  The  Distributed  Monte  Carlo  Fusion  Algorithm 


or 


Ns 

4(4 


1=1 


4-i  (4) 
4-i  (xt) 


n 

jeMi 


(  4-i  41) 
\4-i(4) 


7 

S(x 


xi)(  10) 


From  the  perspective  of  sensor  i  the  Monte  Carlo  dis¬ 
tributed  fusion  algorithm  is: 

Step  0:  Initialisation 

Pick  7  according  to  0  <  7  <  l/max({|A/)|  :  i  G  V }),  and 
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select  a  stopping  time,  fcstop-  For  £  =  1, ...  ,NS  sample 
xkt  ~  nb’  where  irl0  is  a  given  initial  probability  density; 

e-g-  4=0  oc  p(x)g(y.i\x)n  or  7 r|=0  a  p(x)g(yi\x). 

Step  1:  Send,  Receive  ( k  >  0) 

For  each  j  £  Mi  send  the  unweighted  samples  set 
{xk/}^=i  to  agent  j.  For  each  j  £  Mi  receive  {xJk 
Step  2:  Calculate  weights 

If  k  <  kstop  then  for  t  =  1 . . .  Ns  calculate  the  fused 
particle  weight  wlk+1  g  via  (15).  If  k  =  ks top  then  for 

£=1...NS  assign  4+1/  =  4(4 ,4 

Step  3:  Resample 

Resample  with  replacement  Ns  points  xlk+1 1  from  the 

empirical  density  4+1  (4  in  (14).  This  is  equivalent 
to  sampling  from  a  multinomial  distribution  defined  by 
xlk  £  and  w‘k+1 1  If  k  =  ks top  then  the  resampled  set 

{xk+i  4£i  defines  an  unweighted  empirical  probability 
density  of  the  form  (17)  which  is  the  Monte  Carlo 
approximation  of  either  Q  oc  p[x)  IX, ev  g(Vi \x)1/n  or 
Q'  oc  p{x)  nigy  g(yi\x)  (depending  on  the  definition  of 
the  initial  7Tq). 

Note  that  we  have  restricted  the  number  of  iterations  by 
the  parameter  fcstop-  This  is  a  practical  approximation 
obviously  and  one  could  even  use  some  defined  threshold 
on  the  weights  such  that  the  weights  in  (15)  approach 
one  the  algorithm  halts.  Again,  depending  on  whether  the 
initial  likelihoods  are  conditionally  independent  or  not,  one 
may  opt  for  convergence  to  Q  oc  p{x)\\ieV  g{y.i\x)1'n  or 
Q'  <xp(x)Y\zeVg{yi\x). 

2.3  Illustrative  Examples 

In  this  subsection  we  highlight  the  performance  of  the 
distributed  Monte  Carlo  fusion  algorithm  initialised  sim¬ 
ply  by  4=0  oc  p(x)g(yi\x)n  and  where  we  seek  conver¬ 
gence  to  the  true  optimal  Bayesian  fusion  result  Q'  oc 

p(x)  XXgV  9{Vi\x)- 

Firstly,  we  consider  a  random  network  with  5  agents 
shown  in  Figure  1.  The  true  initial  density  4=0  at  each 
agent  was  a  randomly  generated  Gaussian  mixture  with  5 
components.  This  initial  density  was  then  sampled  at  1000 
points  to  generate  each  agent’s  initial  sampled  density 
4=o  (4  which  all  the  agents  would  actually  have  access 
to  in  practice.  The  true  initial  densities  at  each  agent  and 
the  corresponding  random  samples  are  shown  in  Figure  2. 


Fig.  1.  The  network  topology  with  n  =  5  agents. 

The  centralized  Bayesian  solution  (which  is  computable 
from  the  true  initial  Gaussian  mixtures)  was  computed  for 
comparison  and  is  shown  in  Figure  3.  Note  this  cannot  be 
computed  in  practice  because  the  true  initial  (continuous) 
densities  are  unknown  (only  the  samples  are  known)  and 


Fig.  2.  The  actual  initial  densities  are  Gaussian  mixtures 
and  the  samples  are  shown  above  these  for  clarity. 

of  course  we  are  interested  in  distributed  computation 
in  this  work.  However,  this  centralised  solution  is  the 
ideal  solution  and  represents  the  outcome  we  seek  through 
iterative  means  and  via  Monte  Carlo  approximation. 


Fig.  3.  The  ideal  centralised  optimal  Bayesian  fusion  result 
computed  on  the  true  underlying  initial  continuous 
densities.  This  is  not  computable  in  practice  as  we 
suppose  only  the  empirical  (sampled)  density  is  ini¬ 
tially  known  (and  also  the  underlying  true  distribu¬ 
tions  are  unlikely  to  Gaussian  mixtures). 

The  centralised  optimal  Bayesian  solution  was  compared 
to  the  distributed  Monte  Carlo  algorithm  for  Bayesian 
fusion.  The  Monte  Carlo  solutions  are  shown  (in  colour) 
against  the  centralised  solution  (in  black)  in  Figure  4. 


Fig.  4.  The  outcome  from  the  distributed  Monte  Carlo 
fusion  algorithm  after  50  iterations.  The  empirical 
measure  sample  points  are  shown  along  with  continu¬ 
ous  Kernel  density  estimates  generated  from  such  (for 
visualisation  only).  The  centralised  optimal  Bayesian 
solution  shown  in  Figure  3  is  also  shown  again  here 
and  we  can  see  the  Kernel  estimates  of  the  empirical 
(sampled)  densities  closely  approximate  the  optimal 
centralised  solution.  They  also  have  converged  to¬ 
gether  and  reached  consensus. 

From  Figure  4  we  clearly  see  that  each  agent  has  converged 
to  a  common  value  (i.e.  they  have  reached  consensus) 
and  that  this  common  value  has  converged  to  the  desired 
(centralised  Bayesian)  solution. 

We  consider  now  a  similar  illustration  with  10  agents 
and  a  network  with  cycles  as  shown  in  Figure  5.  The 
true  initial  density  7 fk=0  at  each  agent  was  a  randomly 
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generated  Gaussian  mixture  with  5  components.  This 
initial  density  was  then  sampled  at  1000  points  to  generate 
each  agent’s  initial  sampled  density  nlk=0(x)  which  all 
the  agent  would  actually  have  access  to  in  practice.  The 
true  initial  densities  at  each  agent  and  the  corresponding 
random  samples  are  also  shown  in  Figure  5. 


Fig.  5.  The  network  topology  with  n  =  10  agents  and  the 
actual  initial  densities  are  Gaussian  mixtures  and  the 
samples  are  shown  above  these  for  clarity. 

The  centralized  Bayesian  solution  (which  is  computable 
from  the  true  initial  Gaussian  mixtures)  was  computed  for 
comparison  and  is  shown  in  Figure  6.  Again  this  is  just  for 
comparison  and  is  not  realistically  computable  in  practice. 


Fig.  6.  The  ideal  centralised  optimal  Bayesian  fusion  result 
computed  on  the  true  underlying  initial  continuous 
densities.  This  is  not  typically  computable  in  practice. 

The  centralised  optimal  Bayesian  solution  was  compared 
to  the  distributed  Monte  Carlo  algorithm  for  Bayesian 
fusion.  The  Monte  Carlo  solutions  are  shown  (in  colour) 
against  the  centralised  solution  (in  black)  in  Figure  7. 


Fig.  7.  The  outcome  from  the  distributed  Monte  Carlo 
fusion  algorithm  after  100  iterations.  The  empirical 
measure  sample  points  are  shown  along  with  continu¬ 
ous  Kernel  density  estimates  generated  from  such  (for 
visualisation  only).  The  centralised  optimal  Bayesian 
solution  shown  in  Figure  6  is  also  shown  again  here 
and  we  can  see  the  Kernel  estimates  of  the  empirical 
(sampled)  densities  closely  approximate  the  optimal 
centralised  solution.  They  also  have  converged  to¬ 
gether  and  reached  consensus. 

From  Figure  7  we  clearly  see  that  each  agent  has  converged 
to  a  common  value  (i.e.  they  have  reached  consensus) 


and  that  this  common  value  has  converged  closely  to 
the  desired  (centralised  Bayesian)  solution.  This  multi¬ 
modal  fusion  result  shows  the  accuracy  and  potential  of 
the  distributed  Monte  Carlo  fusion  algorithm  (as  only  the 
initial  sample  points  in  Fig  5  are  used  in  initialising  the 
algorithm  and  no  knowledge  about  the  underlying  initial 
continuous  densities  is  assumed). 

3.  A  NOVEL  METHOD  FOR  DISTRIBUTED 
PARTICLE  FILTERING 

In  this  section  we  apply  the  Monte  Carlo  data  fusion 
algorithm  to  systems  of  networked  particle  filters.  Particle 
filters  are  Bayesian  filters  that  represent  their  posterior 
distributions  by  sets  of  unweighted  samples. 

The  standard  algorithm  for  particle  filtering,  also  known 
as  the  bootstrap  filter,  works  as  follows: 

3.1  The  Bootstrap  Filter  Algorithm 

Step  0:  Initialisation  (t  =  0) 

For  £=  1 ...  Ns  sample  x\  t  ~  p(xo|xo),  where  p(xo|xo) 
is  a  known  initial  particle  density. 

Step  1:  Importance  Sampling  (t  >  0) 

For  £  =  l...Ns,  propagate  the  samples  forward  in 
time  through  the  update  (system)  equation,  xlte  = 
f{xlt_1  g,  u\  f)  where  u\  t  is  a  sample  from  the  distribu¬ 
tion  of  process  noise. 

For  t  =  l...iVs,  evaluate  the  importance  weights 
wlt  e  oc  g(y\\x\)  and  normalise,  wt  e  =  1- 

Step  2:  Resample 

For  t  =  1 ...  Ns,  sample  x\  ^  from  the  multivariate 
distribution  p{x\\x\_x,yl)  that  is  constructed  from  the 
samples  xlt  (  and  the  weights  w\  e. 

After  resampling,  the  samples  {x^}^^  define  the  approx¬ 
imation  of  the  posterior  density  of  agent  i.  That  is,  we  are 
approximating  the  true  posterior  given  by 

p(xt\xt-uyl )  ocp(zt|xt_i)sr(2/J|xt)  (18) 

via  a  normalised  empirical  probability  distribution  defined 
by  the  samples  {xt,e 

3.2  The  Distributed  Particle  Filtering  Algorithm 

Consider  a  group  of  agents  (implementing  individual  par¬ 
ticle  filters)  and  an  information  sharing  network  defined 
by  an  undirected  graph  Q  as  defined  in  Section  2. 

The  agents  are  tasked  with  estimating  a  target  state.  Each 
agent  i  £  V,  can  make  observations  y\  about  the  target  at 
time  t. 

Our  goal  in  distributed  particle  filtering  is  to  exchange 
and  combine  posteriors  such  that  each  agent  arrives  at  an 
approximation  of  the  centralised  posterior, 

p(xt\xt-i,ylt)  cxp(xt\xt-i)  Y[g(yl\xt),  (19) 

iev 

provided  that  the  likelihood  functions  g(ylt \xt)  for  all  i  G  V 
are  conditionally  independent. 

To  achieve  this  goal  we  propose  that  the  Monte  Carlo 
fusion  algorithm  is  run  after  Step  2  of  the  bootstrap 
algorithm. 
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From  the  perspective  of  sensor  i: 

Step  0:  Initialisation  (t  =  0) 

Pick  7  according  to  0  <  7  <  l/max({|jV}|  :  j  e  V}). 
For  i  =  1 . . .  Ns  sample  x\t  ~  p(x'0 |cCq),  where  p(x l0\xl0) 
is  a  known  initial  particle  density. 

Step  1:  Importance  Sampling  ( t  >  0) 

For  £  =  1 ...  Ns,  propagate  the  samples  forward  in  time 
x\  t  =  f{x\_lt,u\f)  where  u\  e  is  a  sample  from  the 
distribution  of  process  noise. 

For  £  =  1  ...Ns,  evaluate  the  importance  weights 
wlte  oc  g{y\\x\)'v\  and  normalise. 

Step  2:  Resample 

For  £  =  l...Ns,  sample  x\  e  from  the  multinomial 
distribution  p  constructed  from  the  samples  {xlt 
and  the  weights  {wzt  f}^lv 
Step  3:  Distributed  Ftision 

Run  fusion  algorithm  of  Section  2.2  for  k  =  0 . . .  fcstop. 
For  l  =  1 . . .  Ns  assign  x\^  =  x and  w\  t  =  w{etop>e. 
Note  that  the  weights  should  be  equal  to  1/NS  as  the 
samples  xlk  ^  should  be  unweighted. 


3.3  Illustrative  Examples 


In  this  subsection  we  demonstrate  the  distributed  particle 
filtering  algorithm  on  the  benchmark  nonlinear  dynamical 
system  of  Gordon  et  al.  (1993).  The  process  equation  of 
this  system  is 

xt  =  0.5a:t_i  +  25x*~ 1  +  8cos(1.2(i  -  1))  +  ut,  (20) 
1  +  xt- 1- 


where  Ut  is  Gaussian  white  noise  with  variance,  a2  =  10; 
and  the  sensor  equation  is 


Vt  = 


xt 

20 


(21) 


where  v\  is  Gaussian  white  noise  with  a  random  variance. 
We  initialised  the  filters  with  the  state  xl0  =  0.1  for  all 
i  €  V.  This  represents  a  system  with  a  known  initial  state. 


The  first  simulation  was  conducted  over  a  network  of  10 
agents,  see  Figure  8,  with  Ns  =  1500  samples.  The  true 
state  of  the  system  is  also  shown  in  Figure  8.  We  use 
the  version  of  the  fusion  algorithm  that  converges  to  Q' 
(see  Corollary  2)  as  the  likelihood  functions  are  known 
to  be  conditionally  independent.  Therefore  we  anticipate 
that  the  fusion  algorithm  will  converge  to  the  centralised 
posterior  density. 


Fig.  8.  The  network  topology  of  the  second  distributed 
particle  filter  simulation  with  n  =  10  agents  and  the 
true  state  of  the  nonlinear  dynamical  system  which 
the  agents  are  trying  to  estimate. 


Between  measurements  each  agent  runs  the  distributed 
fusion  algorithm  with  kstop  =  15.  In  Figure  9  we  compare 
the  resulting  estimates  (which  we  call  the  fusion  estimates) 
to  the  desired  state  and  the  centralised  estimate  (com¬ 
puted  at  a  hypothetical  agent  with  access  to  each  agent’s 
likelihood  functions).  The  fusion  estimates  are  shown  by 
the  coloured  dotted  lines  while  the  centralised  estimate 
is  shown  in  solid  blue  and  the  true  state  in  solid  black. 
It  is  clear  that  after  fusion  the  agents  arrive  at  estimates 
that  are  close  to  the  centralised  solution.  In  Figure  10  we 
show  the  estimates  that  would  be  obtained  if  each  node 
ran  the  bootstrap  filter  algorithm  in  isolation  (i.e.  using 
only  their  individual  local  likelihood  functions).  We  refer 
to  these  as  the  isolated  estimates  as  they  do  not  involve  any 
communication  (exchange  of  information).  By  comparison 
to  the  fusion  estimates,  the  isolated  estimates  are  scattered 
around  the  centralised  solution  as  expected. 

In  Figure  11  we  show  the  estimates  generated  by  a  local 
particle  filter  running  at  each  agent  that  makes  use  of  the 
just  the  independent  likelihood  functions  from  each  agent’s 
neighbours.  We  refer  to  these  as  the  local  estimates  and 
this  method.  The  local  estimate  of  i  can  be  thought  of  as 
centralised  estimate  in  the  neighbourhood  of  agent  j.  In  a 
complete  network  the  local  estimates  are  equivalent  to  the 
centralised  estimates. 


Fig.  9.  The  estimates  made  by  the  distributed  particle 
filtering  algorithm  (coloured  dotted  lines).  These  are 
shown  with  the  centralised  estimates  (solid  blue)  and 
the  true  state  (solid  black) .  The  figure  shows  that  the 
agents  are  approximately  reaching  a  consensus  about 
the  fused  posterior  near  the  centralised  solution  as 
desired. 


Fig.  10.  The  isolated  estimates  (coloured  dotted  lines)  are 
shown  along  with  the  centralised  estimate  and  the 
true  state.  As  expected  the  centralised  estimate  falls 
between  the  isolated  estimates  which  are  computed 
by  the  standard  bootstrap  filter  algorithm. 

In  the  second  example  we  randomly  generate  another  net¬ 
work  of  10  agents  and  a  new  underlying  system  trajectory; 


8686 


19th  IFAC  World  Congress 

Cape  Town,  South  Africa.  August  24-29,  2014 


Fig.  11.  The  local  estimates  (coloured  dotted  lines)  are 
shown  along  with  the  centralised  estimate  and  the 
true  state.  The  local  estimates  are  made  by  exchang¬ 
ing  likelihoods  between  neighbours  only.  This  is  im¬ 
portant  as  a  benchmark  for  our  proposed  algorithm 
as  this  is  the  simplest  algorithm  that  may  be  imple¬ 
mented  in  a  network  of  particle  filters. 


True,  Central,  Fusion 


Fig.  13.  The  estimates  made  by  the  distributed  particle 
filtering  algorithm  (coloured  dotted  lines).  These  are 
shown  with  the  centralised  estimates  (solid  blue)  and 
the  true  state  (solid  black) .  The  figure  shows  that  the 
agents  are  approximately  reaching  a  consensus  about 
the  fused  posterior  near  the  centralised  solution  as 
desired. 


see  Figure  12.  Again  the  likelihood  functions  are  condi¬ 
tionally  independent,  however  we  unnecessarily  use  the 
conservative  fusion  method  (for  demonstration  purposes). 
As  the  observations  are  independent  this  method  leads  to 
sub-optimal  estimates. 


Fig.  12.  The  network  topology  of  the  second  distributed 
particle  filter  simulation  with  n  =  10  agents  and  the 
true  state  of  the  nonlinear  dynamical  system  which 
the  agents  are  trying  to  estimate. 

The  agents  run  the  conservative  distributed  fusion  algo¬ 
rithm  with  kstop  =  20  which  converges  towards  the  con¬ 
sensus  value  Q  oc  p{x)  Iljev  g(yj\x)1/n.  This  is  the  con¬ 
servative  posterior.  The  fusion  estimates  were  compared 
to  the  true  state  and  the  centralised  estimate  in  Figure 
13.  We  can  see  that  the  agents  are  approaching  a  common 
value  (i.e.  they  are  close  to  consensus). 

In  Figures  14  and  15  we  show  the  isolated  and  local 
estimates  respectively  (see  the  coloured  dotted  lines).  It 
can  be  seen  that  although  the  local  estimates  are  an 
improvement  on  the  isolated  estimates,  the  conservative 
fusion  estimates  are  better  again.  This  simulation  shows 
that  it  is  possible  to  make  useful  state  estimates  even  when 
the  dependence  between  the  agents  is  unknown  [Bailey 
et  al.  (2012)]. 

4.  CONCLUSION 


Fig.  14.  The  isolated  estimates  (coloured  dotted  lines)  are 
shown  along  with  the  centralised  estimate  and  the 
true  state.  As  expected  the  centralised  estimate  falls 
between  the  isolated  estimates  which  are  computed 
by  the  standard  bootstrap  filter  algorithm. 


Fig.  15.  The  local  estimates  (coloured  dotted  lines)  are 
shown  along  with  the  centralised  estimate  and  the 
true  state.  The  local  estimates  are  made  by  exchang¬ 
ing  likelihoods  between  neighbours  only. 

In  the  future  we  would  like  to  experiment  with  other 
consensus  algorithms  such  as  the  dynamic  consensus  al¬ 
gorithm  of  Zhu  and  Martinez  (2010)  to  reduce  the  com¬ 
munication  overhead. 


We  presented  a  practical  method  for  distributed  data 
fusion  and  particle  filtering.  The  algorithm  is  consensus- 
based  and  involves  the  exchange  of  posteriors  between 
neighbours.  The  proposed  solution  is  robust  to  time- 
varying  network  topologies  including  those  with  cycles  and 
is  guaranteed  consistent. 
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